Universality in the full counting statistics of trapped fermions 
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We study the distribution of particle number in extended subsystems of a one-dimensional non- 
interacting Fermi gas confined in a potential well at zero temperature. Universal features are 
identified in the scaled bulk and edge regions of the trapped gas where the full counting statistics 
are given by the corresponding limits of the eigenvalue statistics in Gaussian unitary random matrix 
ensembles. The universal linnting behavior is confirmed by the bulk and edge scaling of the particle 
number fluctuations and the entanglement entropy. 



The techniques for the manipulation of ultracold 
atomic gases have undergone a rapid development in the 
last decade and have provided experimental access to var- 
ious interesting aspects of many-body quantum systems 
[l|,|2|. The common feature in the experiments is the pres- 
ence of trapping potentials that can be tuned to confine 
particles into effective one-dimensional geometries [3|-|5| 
and theoretical predictions on various properties of ID 
quantum gases p-lsl can directly be tested. 

The strongly correlated phases of ID quantum gases 
are characterized by the simultaneous presence of ther- 
mal and quantum noise. At ultra low temperatures, the 
dominating quantum noise reveals important information 
on the non-local character of correlations in the corre- 
sponding many-body states. In particular, ID Bose liq- 
uids can be probed by measuring the full distribution 
of interference amplitude in experiments [9| , showing re- 
markable agreement with the predictions of the theory 
[lO|. In case of a Fermi gas, an analogous concept is 
the full counting statistics (FCS) [11| which encodes the 
distribution of particle number in extended subsystems. 
The FCS shows interesting properties in the ground state 
of the Fermi gas and, in the non-interacting case, can also 
be used to extract the entanglement entropy of the sub- 
system 12|, |l3[. 



The presence of trapping potentials leaves character- 
istic signatures on the FCS. For the ground state of the 
non-interacting Fermi gas, the FCS was studied in the 
presence of a periodic potential [14| and recently the ef- 
fect of harmonic traps has been analyzed on the particle 
number fluctuations and entanglement [l5|. However, 
the question whether some properties of the FCS hold 
irrespectively of the details of the potential has not yet 
been addressed. 

Here we point out a remarkable universality and show 
that, for a broad class of trapping potentials, the proper 
scaling limits of the FCS in the bulk and edge regime of 
the trapped gas are given by the corresponding eigenvalue 
statistics of Gaussian unitary random matrix ensembles. 
Although the problem is formulated for the free Fermi 
gas, using the Fermi-Bose mapping [l6| the results im- 
mediately carry over to the bosonic Tonks-Girardeaugas 
and thus are relevant for cold-atom experiments [3, 4^ . 



The FCS is defined through the generating function 



X(A) = {g^iXNa)), 



Na^ / p{x)Ax (1) 

Ja 



where Na is the total number of particles in subsystem 
A, given by the integral of the density p{x), and the 
expectation value is taken with the A^-particle ground 
state of the system. For the spinless free Fermi gas, the 
FCS can be expressed as a Fredholm determinant [13J,|l7| 



X(A) = det [1 + (e*^ - 1)Ka\ 
where Ka is an integral operator with the kernel 



N-l 



KA{x,y) = Y^ ipl{x)Lpk{y) 



(2) 



(3) 
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given by the two-point correlation function restricted to 
the domain x,y ^ A. It is constructed from the single- 
particle eigenfunctions of the Schrodinger equation 



2 da;2 



+ [Ek-V{x)]^k{x)^Q 



(4) 



where we we have set fi = m = 1. For simplicity, we con- 
sider a symmetric V{—x) — V{x) and monotonously in- 
creasing trapping potential such that for a; — > cxi one has 
V{x) — >■ oo. The spectrum is thus discrete and for each 
k the solutions ipk{x) admit two classical turning points 
given by the condition V{zLxok) — Ek- In the following 
we will consider x > since the real valued wavefunctions 
must obey the symmetry ipk{—x) = {—l)''(pk{x). 

In the classically allowed region, x < xok, the wave- 
functions are oscillatory with exactly k nodes whereas 
for X > XQk they must vanish exponentially. The wave- 
function which approximates the exact solution on both 
sides is known as the uniform Airy approximation and 
can be derived from a semi-classical treatment of the 
Schrodinger equation [ISJ. Up to the normalization fac- 
tor Cfc, it is given by 



^Pkix) 



Ck 



^/W) 



Ai[±efc(x)] 



(5) 



where the + (— ) sign apphes in the classically forbidden 
(allowed) region and the argument of the Airy function 
is given through 



£.k{x) 



-|2/3 



Pk{z)dz 



(6) 



with xi = niin(a::, xofc) and X2 — niax(a;, xofe). The 
momentum in the integrand of Eq. ^ is defined as 
Pk{x) = y^2\Ek — V{x)\ and the approximate energy 
niveaus Ek can be obtained from the Bohr-Sommerfeld 
quantization formula |19| 



^Ok 



Pk{z)dz = {k + l/2)Tr. 



(7) 



In particular, we will be interested in power-law po- 
tentials of the form Vp{x) = x^ /p with some even integer 
p. Note, that we have set the characteristic length scale 
of the trap to one, which can easily be restored using the 
arguments of trap size scaling [20| . Then the integral in 
Eq. © yields 



Ek^[Mp{k + l/2)] 



- N,^ f r(3/2 + l/,) 

'^ V2pi/P r(i + 1/p) ^ ^ 



with the exponent given hy 6 — p/{p + 2). 

In general, the approximate wavefunctions tpk (x) given 
by Eq. ^ rely on a semi-classical argument and are 
thus expected to reproduce the exact ones ipk{x) only for 
fc ^ 1. In fact, however, the uniform Airy approximation 
gives very good results even for the eigenfunctions of the 
low- lying levels. This is demonstrated on Fig. [ijfor the 
quartic potential V4{x). The eigenfunctions (pk{x) are 
calculated to a high precision by Numerov's method [21 1 
and compared to ipk (x) where the integrals in Eq. ([SJ can 
be given through special (hypergeometric and incomplete 
beta) functions and evaluated numerically. Even though 
the overlap is poor for k = 0, the deviations are already 
very small for fc = 3 and the two functions are essentially 
indistinguishable for fc = 10. 

Far away from the turning point XQk , the uniform Airy 
approximation ([5]) reproduces the WKB-approximation 
|18| which can be obtained from asymptotic expansions 
of Ai[±^fc(x)]. In particular, in the classically allowed 
region one has 



il^kix) 



Ck 



V'^Pkix) 



Pk{z)dz- 



(9) 



For simple potential wells with two classical turning 
points, the normalization constant can be fixed by im- 
posing [22] 



ct 



Pf^ (z)dz = 27r. 



(10) 



Differentiating Eq. ([7]) with respect to fc, one arrives to 
the simple formula C^ = 2^^ and thus the normaliza- 
tion factor accounts for the spectral density. 




FIG. 1. (color online) Exact wavefunctions ipk{x) (symbols) 
and their uniform Airy approximations il^k{x) (lines) for the 
quartic potential Viix). 



The WKB-form of the wavefunctions ^ gives a good 
approximation in the entire classically allowed regime not 
too close to the turning points, but depends on the details 
of the potential V{x). To find universal features of the 
ECS, we first focus on the bulk of the trapped gas and 
we choose the subsystem as the interval A = [—£, (] deep 
in the high density region, (. <C xqn- Expanding the 
functions Pk{x) around x = 0, one obtains 

Ck 



i^kix) 



: cos 



k--^/2Ei:x 



(11) 



which is valid up to linear terms in x. Substituting into 
Eq. ([3]), one arrives to the following sum 



KA(x,y) 




■ V^Ekix - y) 



cos 



(12) 



We are interested in the A^ — > 00 limit of the ECS, 
where the first sum in Eq. (J12p diverges. Therefore we 
introduce new variables and define the scaling limit of 
the kernel as 



Kr{u, v) = lim 



1 



JV^oc ^2F 



--Ka 



N 



\/2Ei\f \/2E 



N 



(13) 



where the subscript refers to the domain [— r, r] of the 
kernel in the scaled variables with the effective length 
r = £^/2En- Introducing the variable z — y/Ek/E^, the 



first sum in (|12p can be converted into an integral while 
the second, alternating sum vanishes in the scaling limit. 
Writing dz = ^^/V^EkEN and using the expression of 
C| in terms of the spectral density, one finds 



Kr{u,v) — — / dzcosz(u — w) 



sin(u — v) 
tt{u — v) 



(14) 



Hence, in the bulk scaling limit, we recover the sine kernel 
which appears in the theory of GUE random matrices. 
Indeed, the probability E{n,t) of finding n eigenvalues 
in an interval [— r, r] in the bulk of the GUE spectrum is 
given by (23| 



E{n,r) 



(-1)" d^ 



dz 



- det(l - zKr) 



(15) 



Fourier transforming Eq. (jl5|) yields the determinant 
in Eq. © with Ka = Kr and thus the bulk FCS of 
the Fermi gas is identical to the bulk GUE eigenvalue 
statistics. 

The bulk scaling limit can be tested through the cu- 
mulants k^ = {—idxy^lnxWlx^Q of the particle num- 
ber. In particular, we calculated the fluctuations K2 as 
a function of £ for the potential V4{x). The numerics is 
simplified by considering, instead of the integral operator 
Ka, the overlap matrix Ca with elements 



Ca. 



kl 



dx(pk{x)(pi{x) 



(16) 



and using TrivT^ = Tr C^ [24|. Then the FCS of Eq. 
^ can be treated as a regular determinant and one has 
K2{£) = TrC^(l — Ca)- This is then compared to ex- 
isting results derived using the asymptotics of Fredholm 
determinants with the sine kernel [17| 

K2{r) = Trif,(l - Kr) = ^(log4r -f 7 + 1) (17) 

where 7 is the Euler constant and we neglected terms 
vanishing for A^ — > 00. 

The results are shown on Fig. [2] for various N with 
the dashed lines representing the curves n,2[r). One can 
see a good agreement for small i but the solid curves 
K2(^) deviate from the scaling prediction as soon as the 
segment size exceeds the size of the flat region in the 
densities p(x) = Ka(x,x), shown on the inset. However, 
the amplitude of the Oix^) corrections to Eq. P^ is 
proportional to y" (0)/-Efe which vanishes for Vp{x) with 
p > 4 and thus the flat density region extends for higher 
p. Note also the strong oscillations in n-iit) as well as 
in p{x) that are results of the alternating sum in Eq. 
([T^ and diminish for higher N . The bulk scaling was 
further tested by calculating the entanglement entropy 
S(t) = -Tr {Ca InC^ + (1 - C^) ln(l - Ca]J and com- 
paring it to the scaling prediction S{r) [2J, |25| with sim- 
ilarly looking results as in Fig. [2l 

The other regime where universal features are expected 
to emerge is near the edge of the high density region. 
Close to the classical turning point, the argument ^fc(x) 
of the Airy function in ([S]) can be expanded around a;ofc 
and yields |18| 



i'kix) 



Ck 



Ai[ak[x - xok)] 



(18) 



K, 0.4 




FIG. 2. (color online) Particle number fluctuations for the 
quartic potential Viix) in the interval [—1, i] (solid lines) com- 
pared to the prediction of Eq. (|17|) in the bulk scaling limit 
(dashed lines) for various A'^. The inset shows the correspond- 
ing density profiles p{x). 



with ak = {2V'{xok))^^^ giving the inverse of the typical 
length scale. The subsystem is now fixed as the interval 
A = {xoN + s/ajyf, 00) starting close to the edge of the 
high-density region and extending to infinity. Note, that 
this choice of the interval A strongly limits the terms con- 
tributing to the sum in Eq. Q since the Airy functions 
in Eq. p^ are shifted gradually to the left for decreasing 
k and for \xok — a^oivl ^ Isj/aAr they become exponen- 
tially small in A. The edge scaling limit of the kernel is 
then defined as 



Ks{u,v) — lim Ka 



Xqn 



OiN 



,XoN 



— (19) 
aw/ 



where the subscript refers to the domain u,v £ (s, 00) in 
the new variables. The factors ak/a^ appearing in the 
arguments of the Airy functions can be approximated as 



ak ^ V (xqn) I 

~ 1 + ^TTTTT \\^^k 

UN W{Xqn) 



Xon) 



(20) 



and since \xQk —xqm\ ^ \s\/o:n, the second term vanishes 
in the limit N ^>- 00 for any well behaved potential. We 
thus set ttfe = QfAT in evaluating (I19p and introduce z = 
aNixQN - XQk). Using xok = V~'^{Ek), one finds dz « 



'Cl/a^f and consequently 



Ks{u,v) 



dzAi{u + z)Ai{v + z) 



Ai(u)Ai'(i;) - Ai'(u)Ai(z;) 



(21) 



Thus we recover the Airy kernel in the edge scaling limit. 
The FCS is then identical to the GUE edge eigenvalue 
statistics [26|, which follows immediately from a formula 
analogous to (jl5|) by replacing r with s. 




N=500 

N=100 

N=50 




FIG. 3. (color online) Rescaled density profiles (dashed lines) 
for V4,{x) near the edge of the high-density region for various 
A''. The uppermost (solid) line shows the N —^ oo scaling 
function p{u), see text. 



FIG. 4. (color online) Particle number fluctuations K2 and 
entropy S in the scaled edge region for various A'^. The solid 
lines show the respective scaling functions K2(s) and S{s). 



As a first test of the scaling limit, we calculated the 
edge density profiles p{x) in the potential V^^x) for vari- 
ous N and compared them to the density scaling function 
p{u) = Ks{u,u) = (Ai'(M))^ - uAi^(u), with the result 
shown on Fig. [31 To reach larger particle numbers N, we 
used, instead of the exact wavefunctions (pk{x), the ipkix) 
in Eq. ([S]), that gives excellent results for the profile in 
the edge region. The finite-size scaling of the data can be 
inferred from the correction term in Eq. (1201) . For power- 
law potentials one has Vp{xoN)/Vp{xoN) = {p — 1)/xqn 
and multiplying by (a^ofc— a;oAr) ~ Oi^ , we obtain a N~^'^ 
scaling of the finite-size corrections, which is consistent 
with the data in Fig. [31 Note, that the exponent of N 
is independent of p, however, the prefactor p — 1 implies 
that the scaling collapse gets worse for larger p which we 
indeed observed in the numerics for p = 6. 

The edge limit of the FCS is verified by calculating 
the particle number fluctuations K2 and entanglement en- 
tropy S of the interval (s, c») in the scaled coordinates. 
They are shown in Fig. [31 for various N and compared to 
the scaling prediction, calculated using a powerful numer- 
ical toolbox for the evaluation of Fredholm determinants 
arising in random matrix theory [27,]. Both K2 and S 
converge slowly to their respective scaling functions. 

In conclusion, we have shown the universality of the 
FCS for a noninteracting trapped Fermi gas in the bulk 
and edge regions, given by the respective eigenvalue 
statistics of GUE random matrices. Interestingly, the 
same universal limits emerge for non- Gaussian unitary 
random matrix ensembles, where the potential V{x) ap- 
pears in the exponential weight function and the semi- 
classical asymptotics of the corresponding orthogonal 
polynomials has to be analyzed 28[ [29[ . This leads to re- 
sults that are surprisingly similar to Eq. ([S]) even though 
the two problems coincide only in the trivial Gaussian 



case, V{x) being the harmonic potential. 

One expects that the same universality would emerge 
for trapped fermions on a lattice and might even gener- 
alize to other potentials. In fact, the connection between 
FCS and GUE statistics has recently been pointed out 
for the time evolution of lattice fermions from a step- 
initial condition [30|. However, in this case the correla- 
tion matrix is unitarily equivalent to the one descri bing 
the ground state of a chain with a linear potential [31 1 
and hence the universality of the FCS carries over to the 
gradient problem. Note that without the lattice, the Airy 
functions in (J18p are the exact eigenstates and the spec- 
trum is continuous, thus the edge limit ([2T|) describes 
the FCS in the entire high-density region while the bulk 
regime is missing. 

It would also be interesting to consider the dynam- 
ical FCS after releasing the gas from the trap. In 
the harmonic case, the corresponding time-dependent 
Schrodinger equation can be solved exactly and the ker- 
nel has, up to an irrelevant phase factor, the equilibrium 
form in appropriately rescaled variables [32|. Thus, the 
limiting scaling forms of the FCS are also unchanged. 
However, for general V{x) the situation is more compli- 
cated and requires a careful analysis. 
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